using System;
using System.Drawing;
using System.IO;
using System.Runtime.InteropServices;

namespace MetaViewAnisotropicDiffusion
{
	namespace Matrix
	{
		public class DoubleMatrix : MetaViewAnisotropicDiffusion.Matrix.MyMatrix
		{
			public double[,] data;
			public DoubleMatrix(int xdim, int ydim)
			{
				this.xDim = xdim;
				this.yDim = ydim;
				this.data = new double[xdim,ydim];
			}
			public DoubleMatrix(int xdim, int ydim, double init)
			{
				this.xDim = xdim;
				this.yDim = ydim;
				this.data = new double[xdim,ydim];
				this.Initialize(init);
			}


			public void Initialize(double init)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] = init;
					}
				}
			}



			public DoubleMatrix Clone()
			{
				DoubleMatrix clone = new DoubleMatrix(this.XDim, this.YDim);
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						clone.data[i,j] = this.data[i,j];
					}
				}
				return clone;
			}

			public void CopyData(DoubleMatrix output)
			{
				// if the sizes doesn't match throw an exception
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						output.data[i,j] = this.data[i,j];
					}
				}

			}


			public void CopyData(Rect inputRect, DoubleMatrix output, int outputX, int outputY)
			{
				// check the boundaries of matrices..
				// Check if the rectangles match..
				// Throw exception
				for(int i=inputRect.startX;i<=inputRect.endX;i++)
				{
					for(int j=inputRect.startY;j<=inputRect.endY;j++)
					{
						int targetX = i + (outputX-inputRect.startX);
						int targetY = j + (outputY-inputRect.startY);

						output.data[targetX,targetY] = this.data[i,j];
					}
				}

			}





			public double Max()
			{
				double max = this.data[0,0];
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						max = this.data[i,j]>max ? this.data[i,j] : max;
					}
				}
				return max;
			}

			public double Min()
			{
				double min = this.data[0,0];
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						min = this.data[i,j]<min ? this.data[i,j] : min;
					}
				}
				return min;
			}


			public double AbsMax()
			{
				double max = Math.Abs(this.data[0,0]);
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						double val = Math.Abs(this.data[i,j]);
						max = val>max ? val : max;
					}
				}
				return max;
			}


			public double AbsMin()
			{
				double min = Math.Abs(this.data[0,0]);
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						double val = Math.Abs(this.data[i,j]);
						min = val<min ? val : min;
					}
				}
				return min;
			}


			public void Abs()
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] = Math.Abs(this.data[i,j]);
					}
				}
			}

			public IntMatrix Sign()
			{
				IntMatrix tmp = new IntMatrix(this.XDim, this.YDim);
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						tmp.data[i,j] = Math.Sign(this.data[i,j]);
					}
				}
				return tmp;
			}

			public double Average()
			{
				double average = 0.0;
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						average += this.data[i,j];
					}
				}
				average /= this.XDim*this.YDim;
				return average;
			}

			public double Sum()
			{
				double average = 0.0;
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						average += this.data[i,j];
					}
				}
				return average;
			}

			public double Std()
			{
				double average = this.Average();
				double std = 0.0;
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						std += (this.data[i,j]-average)*(this.data[i,j]-average);
					}
				}
				std = Math.Sqrt(std/(this.XDim*this.YDim));
				return std;
			}


			public void Sqrt()
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] = Math.Sqrt(this.data[i,j]);
					}
				}
			}

			/// <summary>
			/// Arithmetic operations on matrix elements.
			/// </summary>

			public static DoubleMatrix Add(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				// Check if the dimensions fit together.
				DoubleMatrix output = new DoubleMatrix(mat1.XDim, mat1.YDim);
				for(int i=0;i<mat1.XDim;i++)
				{
					for(int j=0;j<mat1.YDim;j++)
					{
						output.data[i,j] = mat1.data[i,j] + mat2.data[i,j];
					}
				}
				return output;
			}
			public static DoubleMatrix operator+(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				return DoubleMatrix.Add(mat1,mat2);
			}

			public void AddValue(double val)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] += val;
					}
				}
			}

			public static DoubleMatrix operator+(DoubleMatrix mat, double val)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.AddValue(val);
				return tmp;
			}
			public static DoubleMatrix operator+( double val, DoubleMatrix mat)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.AddValue(val);
				return tmp;
			}

			public void AddMatrix(DoubleMatrix mat)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] += mat.data[i,j];
					}
				}
			}

			public static DoubleMatrix Subtract(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				// Check if the dimensions fit together.
				DoubleMatrix output = new DoubleMatrix(mat1.XDim, mat1.YDim);
				for(int i=0;i<mat1.XDim;i++)
				{
					for(int j=0;j<mat1.YDim;j++)
					{
						output.data[i,j] = mat1.data[i,j] - mat2.data[i,j];
					}
				}
				return output;
			}

			public static DoubleMatrix operator-(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				return DoubleMatrix.Subtract(mat1,mat2);
			}
			

			public void SubtractValue(double val)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] -= val;
					}
				}
			}

			public static DoubleMatrix operator-(DoubleMatrix mat, double val)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.SubtractValue(val);
				return tmp;
			}
			public static DoubleMatrix operator-( double val, DoubleMatrix mat)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.SubtractValue(val);
				return tmp;
			}

			public void SubtractMatrix(DoubleMatrix mat)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] -= mat.data[i,j];
					}
				}
			}


			public static DoubleMatrix Multiply(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				// Check if the dimensions fit together.
				DoubleMatrix output = new DoubleMatrix(mat1.XDim, mat1.YDim);
				for(int i=0;i<mat1.XDim;i++)
				{
					for(int j=0;j<mat1.YDim;j++)
					{
						output.data[i,j] = mat1.data[i,j] * mat2.data[i,j];
					}
				}
				return output;
			}

			public static DoubleMatrix operator*(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				return DoubleMatrix.Multiply(mat1,mat2);
			}

			public void MultiplyValue(double val)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] *= val;
					}
				}
			}

			public static DoubleMatrix operator*(DoubleMatrix mat, double val)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.MultiplyValue(val);
				return tmp;
			}
			public static DoubleMatrix operator*( double val, DoubleMatrix mat)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.MultiplyValue(val);
				return tmp;
			}


			public void MultiplyMatrix(DoubleMatrix mat)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] *= mat.data[i,j];
					}
				}
			}


			public static DoubleMatrix Divide(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				// Check if the dimensions fit together.
				DoubleMatrix output = new DoubleMatrix(mat1.XDim, mat1.YDim);
				for(int i=0;i<mat1.XDim;i++)
				{
					for(int j=0;j<mat1.YDim;j++)
					{
						output.data[i,j] = mat1.data[i,j] / mat2.data[i,j];
					}
				}
				return output;
			}

			public static DoubleMatrix operator/(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				return DoubleMatrix.Divide(mat1,mat2);
			}

			public void DivideValue(double val)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						this.data[i,j] /= val;
					}
				}
			}

			public static DoubleMatrix operator/(DoubleMatrix mat, double val)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.DivideValue(val);
				return tmp;
			}
			public static DoubleMatrix operator/( double val, DoubleMatrix mat)
			{
				DoubleMatrix tmp = mat.Clone();
				tmp.DivideValue(val);
				return tmp;
			}


			public void DivideMatrix(DoubleMatrix mat)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						if(mat.data[i,j] != 0d)
						{
							this.data[i,j] /= mat.data[i,j];
						}
					}
				}
			}
			

			public static DoubleMatrix VectorialAdd(DoubleMatrix mat1, DoubleMatrix mat2)
			{
				// Check if the dimensions fit together.
				DoubleMatrix output = new DoubleMatrix(mat1.XDim, mat1.YDim);
				for(int i=0;i<mat1.XDim;i++)
				{
					for(int j=0;j<mat1.YDim;j++)
					{
						output.data[i,j] = Math.Sqrt(mat1.data[i,j]*mat1.data[i,j] + mat2.data[i,j]*mat2.data[i,j]);
					}
				}
				return output;
			}

			public static DoubleMatrix VectorialAdd(ComplexDoubleMatrix mat)
			{
				DoubleMatrix output = new DoubleMatrix(mat.XDim, mat.YDim);
				for(int i=0;i<mat.XDim;i++)
				{
					for(int j=0;j<mat.YDim;j++)
					{
						output.data[i,j] = Math.Sqrt(mat.Real.data[i,j]*mat.Real.data[i,j] + mat.Imag.data[i,j]*mat.Imag.data[i,j]);
					}
				}
				return output;
			}

			public DoubleMatrix SubSample(int ratioX, int ratioY, bool average)
			{

				int nXDim = (this.XDim-1)/ratioX+1;
				int nYDim = (this.YDim-1)/ratioY+1;
				DoubleMatrix tmp = new DoubleMatrix(nXDim, nYDim, 0.0);

				for(int i=0;i<nXDim;i++)
				{
					for(int j=0;j<nYDim;j++)
					{
						if(average)
						{
							double avg = 0.0;
							int count = 0;
							for(int x=0;x<ratioX;x++){
								for(int y=0;y<ratioY;y++)
								{
									if(ratioX*i+x < this.XDim && ratioY*j+y < this.YDim)
									{
										avg += this.data[ratioX*i+x, ratioY*j+y];
										count++;
									}
								}
							}
							tmp.data[i,j] = avg/count;
						}
						else
						{
							tmp.data[i,j] = this.data[ratioX*i, ratioY*j];
						}
					}
				}
				return tmp;

			}

			public void UnitScale()
			{
				double max = this.Max();
				double min = this.Min();
				//init matrix to zero if it is a constant matrix.
				if(max == min)
				{
					this.Initialize(0.0);
				}
				else
				{
					this.SubtractValue(min);
					this.DivideValue(max-min);
				}
			}
			
			public void MeanScale()
			{
				this.DivideValue(this.Average());
			}

			public DoubleMatrix MatrixShift(int shiftX, int shiftY)
			{
				DoubleMatrix shifted = new DoubleMatrix(this.XDim, this.YDim);

				for (int i=shiftX;i<this.XDim;i++)
					for (int j=shiftY;j<this.YDim;j++)
						shifted.data[i-shiftX,j-shiftY] = this.data[i,j];
		
				for (int i=0;i<shiftX;i++)
					for (int j=shiftY;j<this.YDim;j++)
						shifted.data[this.XDim-shiftX+i,j-shiftY] = this.data[i,j];

				for (int i=shiftX;i<this.XDim;i++)
					for (int j=0;j<shiftY;j++)
						shifted.data[i-shiftX,this.YDim-shiftY+j] = this.data[i,j];

				for (int i=0;i<shiftX;i++)
					for (int j=0;j<shiftY;j++)
						shifted.data[this.XDim-shiftX+i,this.YDim-shiftY+j] = this.data[i,j];
			
				return shifted;
			}


			public void ReadRaw(string filename)
			{
				FileStream strm;
				BinaryReader breader;
				try
				{
					strm = new FileStream(filename, FileMode.Open, FileAccess.Read);
					breader = new BinaryReader(strm);
				}
				catch(Exception e)
				{
					Console.WriteLine(e.Message);
					Console.WriteLine("Cannot open" + filename + "for reading");
					return;
				}
				for(int j=0; j<this.yDim; j++)
				{
					for(int i=0; i<this.xDim; i++)
					{
						this.data[i,j] = breader.ReadDouble();
					}
				}
				breader.Close();
				strm.Close();
			}

			public void ReadRawFloat(string filename)
			{
				FileStream strm;
				BinaryReader breader;

				try
				{
					strm = new FileStream(filename, FileMode.Open, FileAccess.Read, FileShare.Read, this.xDim*this.yDim*4);
					breader = new BinaryReader(strm);
				}
				catch(Exception e)
				{
					Console.WriteLine(e.Message);
					Console.WriteLine("Cannot open" + filename + "for reading");
					return;
				}

				
				for(int j=0; j<this.yDim; j++)
				{
					for(int i=0; i<this.xDim; i++)
					{
						this.data[i,j] = (double)breader.ReadSingle();
					}
				}
				breader.Close();
				strm.Close();
			}

			public void WriteRaw(string filename)
			{
				FileStream strm;
				BinaryWriter bwriter;
				try
				{
					strm = new FileStream(filename, FileMode.Create, FileAccess.Write);
					bwriter = new BinaryWriter(strm);
				}
				catch(Exception e)
				{
					Console.WriteLine(e.Message);
					Console.WriteLine("Cannot open" + filename + "for writing");
					return;
				}
				for(int j=0; j<this.yDim; j++)
				{
					for(int i=0; i<this.xDim; i++)
					{
						bwriter.Write(this.data[i,j]);
					}
				}
				bwriter.Close();
				strm.Close();
			}

			public static DoubleMatrix ReadImage(string filename)
			{
				Bitmap image = new Bitmap(filename);
				DoubleMatrix temp = new DoubleMatrix(image.Width, image.Height);
				for(int i=0;i<temp.XDim;i++)
				{
					for(int j=0;j<temp.YDim;j++)
					{
						temp.data[i,j] = 255*(image.GetPixel(i,j)).GetBrightness();
					}
				}
				return temp;
			
			}

			public static DoubleMatrix ReadBitmap(Bitmap image)
			{
				DoubleMatrix temp = new DoubleMatrix(image.Width, image.Height);
				for(int i=0;i<temp.XDim;i++)
				{
					for(int j=0;j<temp.YDim;j++)
					{
						temp.data[i,j] = 255*(image.GetPixel(i,j)).GetBrightness();
					}
				}
				return temp;
			
			}

			public void WriteImage(string filename)
			{
				Bitmap image = new Bitmap(this.XDim, this.YDim);

				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						image.SetPixel(i, j, Color.FromArgb((int)Math.Round(this.data[i,j]), (int)Math.Round(this.data[i,j]), (int)Math.Round(this.data[i,j])));
					}
				}
				
				image.Save(filename);

			}

			public Bitmap WriteBitmap()
			{
				Bitmap image = new Bitmap(this.XDim, this.YDim);

				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						image.SetPixel(i, j, Color.FromArgb((int)Math.Round(this.data[i,j]), (int)Math.Round(this.data[i,j]), (int)Math.Round(this.data[i,j])));
					}
				}
				
				return image;

			}

			public void WriteImageScaled(string filename)
			{
				DoubleMatrix tmp = this.Clone();
				tmp.UnitScale();
				tmp.MultiplyValue(255.0);
				tmp.WriteImage(filename);
			}

			public Bitmap WriteBitmapScaled()
			{
				DoubleMatrix tmp = this.Clone();
				tmp.UnitScale();
				tmp.MultiplyValue(255.0);
				return tmp.WriteBitmap();
			}

			public void Threshold(double th)
			{
				for(int i=0;i<this.XDim;i++)
				{
					for(int j=0;j<this.YDim;j++)
					{
						if(this.data[i,j] >th)
						{
							this.data[i,j] = 1.0;
						}
						else
						{
							this.data[i,j] = 0.0;
						}
					}
				}
			}

			public static explicit operator IntMatrix(DoubleMatrix dMat)
			{
				IntMatrix result = new IntMatrix(dMat.XDim, dMat.YDim);
				for(int i=0;i<dMat.XDim;i++)
				{
					for(int j=0;j<dMat.YDim;j++)
					{
						result.data[i,j] = (int)dMat.data[i,j];
					}
				}
				return result;

			}
			
			public static implicit operator MultiValuedDoubleMatrix(DoubleMatrix dMat)
			{
				MultiValuedDoubleMatrix result = new MultiValuedDoubleMatrix(dMat.XDim, dMat.YDim, 1);
				for(int i=0;i<dMat.XDim;i++)
				{
					for(int j=0;j<dMat.YDim;j++)
					{
						result.comp[0].data[i,j] = dMat.data[i,j];
					}
				}
				return result;

			}


			public static double Dist1(DoubleMatrix m1, DoubleMatrix m2)
			{
				double dist = 0d;
				for(int i=0;i<m1.XDim;i++)
				{
					for(int j=0;j<m1.YDim;j++)
					{
						dist += Math.Abs(m1.data[i,j] - m2.data[i,j]);
					}
				}
				dist /= m1.XDim*m1.YDim;
				return dist;
			}


		}

	}
}









